source("scripts/recursos.R")
source("scripts/funciones.R")4 CVM ppp
4.1 Pretratamientos de la Información
Cargar Recursos Externos
Definición de variables generales
## geographic projections
crs_latlon <- "+proj=longlat +datum=WGS84 +no_defs"
crs_utm <-
"+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"Definición de variables específicas
# List of crimes/delitos and cities/ciudades
crimes = c(
"Minor property offenses",
"Robberies",
"Burglaries",
"Drunkennes, damages and disorders",
"Family violence and aggression",
"Injuries, drugs and weapons",
"High violence and murder"
)
delitos = c("hurto", "robo", "asalto", "amen", "abus", "crim", "viol")
#ciudad de ejemplo
cities <- c("Coquimbo-Serena")
ciudades <- c("urb_4_18")
obs <- 1
delito <- delitos[1]Lectura de Datos Espaciales
Para el caso del presente ejemplo práctico se utiizará la ciudad de Coquimbo-Serena.
city <- ciudades[obs]
cityp <- readRDS(paste0("data/",city,".rds"))
cityp <- spTransform(cityp, CRS(crs_utm))El objeto espacial cityp cuyos puntos corresponde donde ha ucurrido eventos deictuales, cuyo formato SpatialPointDataFrame.
Matriz de vecindad espacial: Se utilizó la función propia spatial_weights() @(sec-spatial-weights)
# matriz de vecindad espacial
nb <- spatial_weights(city_df = cityp, nvec= 12)Ventana de Trabajo para ppp (Point Pattern Plane)
Crea un objeto de clase “ppp” que representa un conjunto de datos de patrones de puntos en el plano bidimensional. Por lo anterior, se debe crear una ventana de observación w.
#make windows
ext <- raster::extent(cityp)
x_min <- ext[1] + 100
x_max <- ext[2] - 100
y_min <- ext[3] + 100
y_max <- ext[4] - 100
w <- as.owin(c(x_min,x_max, y_min, y_max))Transformaciones a los datos.
Se utilizó funciónes de transformacion general mencionadas en Section 6.2 que corresponden a las siguientes:
- Transformación Cúbica
- Transformación Logarítmica
- Cáculo del Lag Espacial
cityp <- cityp %>%
st_as_sf(crs = 4326) %>%
mutate(across(.cols = all_of(delitos),
.fns = t_cub, .names = "cub{.col}")) %>%
mutate(across(.cols = all_of(delitos),
.fns = t_log1p, .names = "log{.col}")) %>%
mutate(across(.cols = all_of(delitos),
.fns = lag_spatial, nb, .names = "lag{.col}")) %>%
as("Spatial")| zona | id | x | y | pflot | viol | crim | abus | amen | robo | hurto | asalto | denspob | educ | tamhog | constr | ofcom | vulzon | mobil | lpflot | neduc | ldenspob | lofcom | lconstr | cubhurto | cubrobo | cubasalto | cubamen | cubabus | cubcrim | cubviol | loghurto | logrobo | logasalto | logamen | logabus | logcrim | logviol | laghurto | lagrobo | lagasalto | lagamen | lagabus | lagcrim | lagviol |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 4101161005 | 22 | -71.239 | -29.862 | 1.000 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 144.756 | 10.867 | 3.170 | 4328.618 | 0.000 | 0.367 | 0.079 | 0.178 | 0.477 | 0.767 | 0.000 | 0.840 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.083 | 0.083 | 0.083 | 0.417 | 0.250 | 0.083 | 0 |
| 4101161005 | 32 | -71.243 | -29.863 | 1.000 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 227.772 | 10.716 | 3.393 | 3830.009 | 0.000 | 0.367 | 0.079 | 0.178 | 0.466 | 0.846 | 0.000 | 0.829 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.000 | 0.083 | 0.083 | 0.917 | 0.583 | 0.083 | 0 |
| 4101161005 | 35 | -71.240 | -29.863 | 1.034 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 171.133 | 11.801 | 3.172 | 4328.618 | 0.000 | 0.367 | 0.079 | 0.182 | 0.550 | 0.796 | 0.000 | 0.840 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0.693 | 0.000 | 0.000 | 0 | 0 | 0.000 | 0.083 | 0.000 | 0.500 | 0.250 | 0.083 | 0 |
| 4101161005 | 36 | -71.239 | -29.863 | 1.000 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 154.603 | 10.870 | 3.135 | 4328.618 | 0.000 | 0.367 | 0.079 | 0.178 | 0.478 | 0.779 | 0.000 | 0.840 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.083 | 0.083 | 0.083 | 0.417 | 0.250 | 0.083 | 0 |
| 4101161005 | 37 | -71.238 | -29.863 | 1.025 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 155.145 | 10.858 | 3.148 | 3230.309 | 0.000 | 0.367 | 0.079 | 0.181 | 0.477 | 0.779 | 0.000 | 0.815 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.083 | 0.083 | 0.083 | 0.000 | 0.250 | 0.083 | 0 |
| 4101161001 | 48 | -71.244 | -29.864 | 1.125 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 168.179 | 11.958 | 3.336 | 2384.222 | 14.682 | 0.292 | 0.122 | 0.193 | 0.562 | 0.793 | 0.574 | 0.788 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.000 | 0.000 | 0.250 | 0.500 | 0.417 | 0.083 | 0 |
| 4101161005 | 49 | -71.243 | -29.864 | 1.071 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 212.743 | 10.978 | 3.415 | 3795.759 | 0.000 | 0.367 | 0.079 | 0.187 | 0.486 | 0.834 | 0.000 | 0.828 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.000 | 0.083 | 0.250 | 0.917 | 0.583 | 0.083 | 0 |
| 4101161005 | 50 | -71.242 | -29.864 | 1.000 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 220.188 | 10.705 | 3.409 | 4026.613 | 0.000 | 0.367 | 0.079 | 0.178 | 0.465 | 0.840 | 0.000 | 0.834 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0.000 | 0.693 | 0.693 | 0 | 0 | 0.000 | 0.083 | 0.000 | 1.000 | 0.417 | 0.083 | 0 |
| 4101161005 | 51 | -71.241 | -29.864 | 1.021 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 243.501 | 10.120 | 3.337 | 4328.618 | 0.000 | 0.367 | 0.079 | 0.180 | 0.419 | 0.857 | 0.000 | 0.840 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0.000 | 0.693 | 0.693 | 0 | 0 | 0.000 | 0.083 | 0.083 | 0.667 | 0.250 | 0.000 | 0 |
| 4101161005 | 52 | -71.240 | -29.864 | 1.022 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 320.023 | 10.720 | 3.324 | 4328.618 | 0.000 | 0.367 | 0.079 | 0.180 | 0.466 | 0.905 | 0.000 | 0.840 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.000 | 0.083 | 0.167 | 0.583 | 0.250 | 0.000 | 0 |
4.2 Evaluación Cramer Von Mises
4.2.1 Modelos MLB
Para transformar los datos espaciales de los delitos y las predicciones de los modelos se utilizó la función propia ppp_maker() Section 6.4 que retorna el objeto ppp.
Lectura del Modelo.
model_mlb <- readRDS(paste0("output/mlb_",city,"_",delito,".rds"))Transformaciones de las predicciones del modelo a ppp.
ppp_mbl_pred <- ppp_maker(model = model_mlb, city = cityp,
delito = delito, ppp_predict = T, w = w)
plot(density(ppp_mbl_pred, adjust=.25))
Transformaciones de los delitos reales a a ppp.
ppp_mbl <- ppp_maker(model = model_mlb, city = cityp,
delito = delito, ppp_predict = F, w= w)
plot(density(ppp_mbl, adjust=.25))
Evaluación cdf.test.()
cvm_mlb <- cdf.test(ppp_mbl, as.im(ppp_mbl), test="cvm")
plot(cvm_mlb)
attr_model <- attr(cvm_mlb, which = "frame")
im <- attr_model$values$Zimage
r_mlb <- raster(im)
proj4string(r_mlb)=crs_utm
values(r_mlb) <- ifelse(values(r_mlb)==0, NA, values(r_mlb))
plot(r_mlb)
pred_mlb <- predict_df(model = model_mlb, city = cityp, longlat = F)
pred_mlb <- st_as_sf(x = pred_mlb,
coords = c("x", "y"),
crs = 32719) %>%
filter(!is.na(delitos))
mapview(r_mlb, na.color= NA)+
mapview(pred_mlb, zcol = "ypred",hide = TRUE, cex = 2)+
mapview(pred_mlb, zcol = "delitos", hide = TRUE, cex = 2)Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
+b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
+wktext +no_defs +type=crs
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
4.2.2 Modelos SAR
city <- ciudades[obs]
model_sar <- readRDS(paste0("output/sar_",city,"_",delito,".rds"))
ppp_sar_pred <- ppp_maker(model = model_sar, city = cityp,
delito = delito, ppp_predict = T,
w = w, nb = nb)This method assumes the response is known - see manual page
plot(density(ppp_sar_pred, adjust=.1))
ppp_sar <- ppp_maker(model = model_sar, city = cityp,
delito = delito, ppp_predict = F,
w = w, nb = nb)
plot(density(ppp_sar, adjust=.25))
# st <- syrjala.test(ppp_mbl, ppp_mbl_pred, nsim = 999)
cvm_sar <- cdf.test(ppp_sar, as.im(ppp_sar_pred), test="cvm")
plot(cvm_sar)
attr_model <- attr(cvm_sar, which = "frame")
im <- attr_model$values$Zimage
r_sar <- raster(im)
proj4string(r_sar)=crs_utm
values(r_sar) <- ifelse(values(r_sar)==0, NA, values(r_sar))
# plot(r_sar)pred_sar <- predict_df(model = model_sar, city = cityp, longlat = F)This method assumes the response is known - see manual page
pred_sar <- st_as_sf(x = pred_sar,
coords = c("x", "y"),
crs = 32719) %>%
filter(!is.na(delitos))
mapview(r_sar, na.color= NA)+
mapview(pred_sar, zcol = "ypred",hide = TRUE, cex = 2)+
mapview(pred_sar, zcol = "delitos", hide = TRUE, cex = 2)